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ABSTRACT 

Context. The formation of resonant pairs of planets in exoplanetary systems involves planetary migration in the protoplanetary disc. 
After a resonant capture, the subsequent migration in this configuration leads to a large increase of planetary eccentricities if no 
damping mechanism is applied. This has led to the conclusion that the migration of resonant planetary systems cannot occur over 
large radial distances and has to be terminated sufficiently rapidly through disc dissipation. 

Aims. In this study, we investigate whether the presence of an inner disc might supply an eccentricity damping of the inner planet, 
and if this effect could explain the observed eccentricities in some systems. 

Methods. To investigate the influence of an inner disc, we first compute hydrodynamic simulations of giant planets orbiting with a 
given eccentricity around an inner gas disc, and measure the effect of the latter on the planetary orbital parameters. We then perform 
detailed long term calculations of the GJ 876 system. We also run N-body simulations with artificial forces on the planets mimicking 
the effects of the inner and outer discs. 

Results. We find that the influence of the inner disc can not be neglected, and that it might be responsible for the observed eccentric- 
ities. In particular, we reproduce quite well the orbital parameters of a few systems engaged in 2: 1 mean motion resonances : GJ 876, 
HD 73 526, HD 82 943 and HD 128 31 1. Finally, we derive analytically the effect that the inner disc should have on the inner planet 
to reach a specific orbital configuration with a given damping effect of the outer disc on the outer planet. 

Conclusions. We conclude that an inner disc, even though difficult to model properly in hydro-dynamical simulations, should be 
taken into account because of its damping effect on the eccentricity of the inner planet. By including this effect, we can explain quite 
naturally the observed orbital elements of the pairs of known resonant exoplanets. 

Key words. Accretion, accretion discs - Planets and satellites : formation - Celestial mechanics 



1. Introduction 

The orbital evolution of a system consisting of very young pro- 
toplanets is governed by disc-planet and mutual gravitational in- 
teractions. In case of differential migration the semi-major axis 
ratio of two planets varies with time and - in the situation of con- 
vergent migration - capture in a resonant configuration may oc- 
cur. Indeed, a large fraction of the observed multi-planet systems 
contain a pair of planets engaged in a resonance. Here, we are in- 
terested in mean motion resonances (MMR) where the ratio of 
the (mean) orbital periods of the outer to the inner planet equals 
that of two small integers. Among the 6 systems known to be in 
a MMR, 4 have a ratio of 2: 1 : GJ 876, HD 82 943, HD 128 31 1 
and HD 73 526. The system discovered first to lie in a 2:1 con- 
figuration (GJ 876) is interesting in several aspects. The planets 
are both very massive (0.56 and 1.94 Mj up ), particularly when 
considering the small mass of the central star of only 0.33 M . 
The relatively short orbital periods of the planets (« 30 and » 60 
days) have allowed for very accurate determination of their or- 
bital elements, which are stated in Table Q] In GJ 876 the two 
outer planets are 'deep' in the 2:1 MMR, i.e. the apsidal lines of 
the two osculating orbital ellipses are always aligned and librate 
with very small amplitudes only (so called apsidal resonance or 
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corotation). As a consequence of the apsidal resonance the plan- 
etary eccentricities show only small variations with time. A res- 
onant configuration like that in GJ 876 can only be established 
through the action of dissipative effects such as disc-planet inter- 
action. In fact, the mere existence of systems engaged in MMRs 
is one of the strongest indication that planetary migration has 
indeed occurred during the early evolution of planetary systems. 

The first detailed modelling of GJ 876 has been conducted by 
iLee & Peald (l2002al) who performed customised 3-body simula- 
tions of a central host star and two planets with additional (dis- 
sipative) forces mimicking the effects of disc-planet interaction. 
In these type of simulations it is assumed that a pair of plan- 
ets is still embedded in the protoplanetary disc, which consists 
of an outer disc only while the inner disc (inside of the inner 
planet) has already been lost through effects like accretion onto 
star and planets or final evaporation. In such a configuration only 
the outer planet is in contact with an even further protoplanetary 
disc and it experiences typically negative Lindblad torques and 
migrates inward. On the contrary the inner planet has no ambi- 
ent material and doe s not feel any disc t orque. In terms of the 3- 
body simulations bv lLee & Peald d2002al) this implies additional 
forces that reduce the semi-major axis and eccentricity of the 
outer planet, while the inner planet feels only the direct gravita- 
tional forces of the star and outer planet. Capture into a resonant 
configuration can now occur when during the inward migration 
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the outer planet crosses the location at which the mean orbital 
periods have a ratio of two small integers. 

The eccentricity damping action of the outer disc onto the 
outer planet is typically parameterised through the migration 
rate, i.e. 



-K 



(1) 



where a and e are the semi-major axis and eccentricity of 
the planet, respectively, while the constant factor K relates the 
damping rates. For low mass planets in the linear regime, it is 
known that the eccentricity damping time scale t e = -eje is 
about h 2 times shorter than the migration time scale t a = -a /a 
(Tanaka & Ward 2004), where h is the aspect ratio of the disc 
(generally a few percent). This short time scale has been con- 
firmed recently through fully n onlinear hydro-dynam ical sim- 
ulations in 2 and 3 dimensions (ICresswell et al.ll2007l) . Hence, 
for low mass planets one may safely assume that they migrate 
inward on nearly circular orbits unless otherwise disturbed by 
additional objects in the system. For higher planetary masses the 
value (and also the sign) of K is not known exactly, but due to 
the opening of a gap it is to be expected that the damping on the 
eccentricity wi ll be reduced with re spect to the linear case. The 
simulations bv lLee & Peald d2002al) indicate now that for values 
of K equal to unity the resonant action leads to a strong growth 
of eccentricity of both the inner and outer planets much larger 
than th e observed values (e\ « 0.03 and e2 ~ 0.22). lLee & Peala 
(2002a) pointed out that the final state of the system is deter- 
mined by K. For a vanishing K, i.e. no eccentricity damping, the 
plane tary eccentricities co ntinue to grow. To match the observa- 
tions [CeejfcPilll d2002al) had to assume a value of K — 100, a 
value that appears to be very large considering the high masses 
of the planets. 

The convergent migration of two massive planets has been 
demonstrated in a variety o f hydro-dynamical simu l ations ( Kiev 
2000t iBrvden et al.1 120001: ISnellgrove etail 1200 U iPapaloizou 
2003). Through multi-dimensional hydro-dynamical simulations 
of resonant planetary systems it has been shown that for masses 
in the Jupiter re gime the value of K lies around 1 - 10 at most 
dKlev et al. l |2004|> . and in a detailed study of the system GJ 876 



Kiev et al. ( 20051) have shown that in hydro-dynamical Simula 
tions where the inner disc has been depleted the final eccentrici- 
ties of the planets will always be much larger than the observed 
values unless one assumes that the outer disc dissipates rapidly 
on the viscous time scale. Hence, this scenario does not allow 
for the migration of the resonant planets over a larger radial 
distance. While in hydro-dynamical simulations of discs with 
embedded planets it is often found that the inner disc is de- 
pleted rapidly, this may in reality not be the case and be an arte- 
fact of inappropriate inner boundary conditions. As shown by 
Crida & MorbideH (120071) . even in the presence of a planet an 
inner disc will survive for much longer than found previously. In 
this situation it can be expected that the inner disc will have a dy- 
namical influence on the inner planet and induces possibly some 
additional damping of the eccentricity. The influence of such 
an eccentricity damping of the inner planet was mentioned first 
bv lLee &Pealel (l2002al) for the particular case of GJ 876. They 
found that a value of K = 10 for both planets yields final ec- 
centricities in the observed range. The first full hydro-dynamical 
study in this direction - including an inner disc - was done for 
the resonant system HD 73 526 bv lSandor et al.l (120071) ; in their 
study, they have shown that the inclusion of an inner disc indeed 
leads to an additional eccentricity damping of the inner planet, 



and allows more extended radial migration with reasonable final 
eccentricities. 

In this paper we analyse this effect in more detail and investi- 
gate the dynamical influence that an inner disc has on a planet. In 
Sect.|2]we perform a sequence of hydrodynamic simulations and 
measure the torque and power exerted by the disc on the planet 
and evaluate the change in eccentricity and semi-major axis as 
a function of the planetary eccentricity. In Sect. [3] we perform 
a full time evolution of a pair of planets embedded in a proto- 
planetary disc for realistic parameters with specific application 
to GJ 876. We show that the torque and the power generated by 
the inner disc yield an effective damping of e which results in 
moderate final eccentricities even for extended radial migration. 
Last, in Sect. [4] we apply an eccentricity damping to the inner 
planet, in the frame of N-body simulations with artificial non- 
conservative forces to mimic the effect of the disc. We apply this 
to a few exoplanetary systems, and try to recover the observed 
orbital elements with a realistic migration scenario. Sect.[5]sum- 
maries our results. 



2. Effect of a gaseous inner disc on the orbital 
elements of a planet on an eccentric orbit 

To measure the influence that an inner disc has on the orbital 
elements of a planet orbiting around it on an eccentric orbit, we 
perform a suite of hydro-dynamical simulations. Here the disc 
is treated as a two-dimensional gas that lies in the orbital plane 
of the planet. The disc material is (basically) only present inside 
the planetary orbit, so that the effect of the inner disc can be 
isolated. The planet has a mass M p — 2.14 X 10~ 3 M» and has a 
fixed orbit with semi major axis a = 1, and a given eccentricity 
e, that changes from one simulation to another. 

The disc is treated as a non-self-gravitating gas that never- 
theless can interact gravitationally with the planet. From the ex- 
pressions for the planetary energy E - -G(M„ + M p )M p /2a and 

angular momentum H - M p ^G(M» + M p )a{\ - e 2 ), one can 
easily derive : 



a/a = —E/E 

e 2 -\lE H 
e/e = —r^r- I — + 2 



2e 2 



H 



(2) 
(3) 



In the present simulations we keep the orbit of the planet fixed 
but monitor the torque (H) and power (E) acting on the planet 
(averaged over one orbit). We have checked that if the planet is 
released from its fixed orbit, its eccentricity and semi major axis 

follow the expected evolution. 

The code we use is FARGO, by iMassetl (l2000allbl) which 
is a 2D hydro-code using cylindrical coordinates (r, if) with an 
isothermal equation of state. Thus, the sound speed is given by 
c s = /irQ, where Q is the local angular velocity, r is the distance 
to the central star, and h is the aspect ratio, which is here 0.05. 
The gas viscosity v is g iven by an a-prescription (v = ac s hr, 
Shak ura & Sunvaevll 19731) . with a = 10~ 2 . The grid covers the 
region from 0.4 to 1 .62 in radius, divided in 1 12 elementary rings 
(logarithmically spaced), themselves divided in 500 sectors. The 
inner boundary condition is non-reflecting. More precisely, at 
every time step the density in the zeroth ring is set equal to the 
one in first ring, rotated by a suitable angle to mimic wave prop- 
agation, which avoids wave reflection ; then, the density in the 
zeroth ring is shifted such that the azimuthally averaged density 
is the same as initially. The outer boundary is open, which means 
that outflow of gas out of the grid is permitted. 
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Fig. 1. Gas density profiles after 400 orbits for a few planet ec- 
centricities, and common initial profile (plain line). 



The coordinate system is centred on the star. It can be non- 
rotating, corotating with the planet, or rotating at a constant an- 
gular velocity with a period equal to that of the planet 



„ / G(M„ + M p ) 
00 = 2 



1/2 



(4) 



Even with the FARGO algorithm, which at each ring transforms 
essentially to a corotating frame, the choice of a different rota- 
tion rate of the coordinates may change the gas dynamics close 
to the planet. We have performed comparison simulations in all 
three frames and display some results below. In the corotating 
frame, gas in the Hill sphere simply rotates around the planet, 
whose motion is only radial from apoastron to periastron and 
back. In the non-rotating frame, the motion of the planet around 
the star in the grid tends to spread artificially the Hill sphere. 
For large eccentricities, however, the corotating frame is not ad- 
vantageous because it is not rotating at a constant angular veloc- 
ity, which implies additional terms in the equations of motion. 
Therefore we prefer the last option, a frame rotating with the 
constant speed Qo in which the planet describes an epicycle. 

The initial density profile corresponds to an approximate gap 
opened by the planet, without the outer disc, so that an equilib- 
rium profile and a stationary regime are quickly reached. The 
initial profile can be seen on Fig. Q] as a solid line. The initial 
total mass of gas in the disk is 1.53 x 10~ 3 M« = 1\%M P . The 
final density profiles after 400 orbits for various planet eccen- 
tricities are also displayed on Fig.Q] The profiles are taken when 
the planet is at apoastron, and the density spike at r — 1 + e cor- 
responds to the Hill sphere of the planet. The larger the eccen- 
tricity, the more depleted the inner disc. The outer disc is clearly 
empty. In the 5 presented cases, the simulation has been run in 
the non-rotating frame. 

2. 1 . Average effect of the inner disc 

The power and the torque from the disc on the planet become al- 
most constant after about 200 orbits in most cases, and we mea- 
sure them after 400 orbits. In the measure of the force exerted 
by the gas on the planet, we exclude material in the Hill sphere, 
using a tapering function given by : 
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Fig. 2. Tapering function f(d) defined by Eq. (O for four values 
of the parameter p. 



where d is the distance to the planet and r# is its Hill radius 
r H = fp(M p /3M,) 1 ^ (r p being the distance between the planet 
and the star) ; p is a dimensionless parameter set equal to 0.8. 
This function is plotted on Fig. [2] for different values of p. We 
stress that the use of such a tapering is really important as the gas 
close to the planets exerts on it a strong force, which should not 
be taken into account because this gas is gravitationally bound 
to the planet and should be considered as a part of it. With self- 
gravity, this gas should feel the same force as the planet and 
naturally follow it ; using a tapering frees the planet from car- 
rying artificially this material. The influence of the shape of the 
tapering function will be discussed below. 

The results are shown in Fig.[3]for the three above mentioned 
options for the rotation rate of the coordinate system. The quan- 
tities are normalised by p — M<n sc /M p , because the force felt by 
the planet is proportional to the gas density. The two top panels 
display the influence of the inner disc on the orbital elements e 
and a. For e > 0.1, the eccentricity is effectively damped on a 
time scale of about two thousands of orbits, as can be seen in 
the bottom left panel displaying T e - —el(elp) ; the dispersion 
in the values of r e for 0.1 < e < 0.15 is simply d ue to the fact 
that e i s clo se to zero. This con firms the idea by iLee & Peald 
(2002b) and Sa ndor et al.l (120071) that an inner disc has indeed a 
damping effect on the planetary eccentricity. For low eccentric- 
ities (e < 0.1), the influence of the disc is small (|e| < 4 x 10~ 5 
orbit 1 ), and could even lead to a small excitation of the eccen- 
tricity. 

It is well-known that a planet o n a circular orbit exerts 
a negative torque on the in ner disc dLin & Papaloizoul Il979t 
iGoldreich & Tremaind fl980). It repels the gas, leading to the 
opening of a gap. When the density gradient at the gap edge is 
steep enough, an equilibrium arises (see for instance TCrida et all 
2006). Then, the planet feels a positive torque from the inner 
disc. For e — 0, we find that the planet feels a positive torque 
(which is equal to the power), as expected ; consequently, a > 0. 
But as e increases, the power decreases, and for e > 0.3, the semi 
major axis variation expected if one releases the planet is nega- 
tive (top right panel). In a sense, the inner disc attracts the planet 
more and more as the eccentricity grows. 

As a consequence of the variations of a and e with e, the 
simple modelling given by Eq. ([TJ with K constant seems to be 
an over simplification. The bottom right panel of Fig. [3]displays 
(e/e)/(d/a) — K as a function of e. It varies a lot. The dispersion 
around e = 0.3 - 0.35 is due to the change of sign of a. For 
0.1 < e < 0.25, a K factor of the order of minus a few seems 
reasonable, but a precise value can not be stated. However, it is 
clear that the effect of the inner disc on the planet can not be 
neglected. In particular the left panels show that a significant 
damping of the eccentricity has to be taken into account. 
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Fig. 3. Influence of the inner disc on a planet on a fixed orbit, 
as a function of the eccentricity. NR (+ symbols) : in the non- 
rotating frame ; CF (x symbols) : in the corotating frame ; QF 
(open squares) : in the frame rotating monotonically with the 
same period as the planet. 



2.2. More detailed analysis 

To investigate the physical mechanism more thoroughly, we fo- 
cus on the case where e = 0.15 and where the frame is rotating 
at a angular constant velocity Q(j. On Fig.|4j we plot the specific 
torque and power felt by the planet during one orbit, starting at 
apoastron. The units are normalised such that a = G = M* = 1, 
and thus Q = Vl.00214. Four curves are presented, corre- 
sponding to measures on the same orbit with four different val- 
ues of the parameter p in Eq. (0, in the range [0.6; 0.9]. The 
larger the p factor, the smaller the amplitude of the oscillation. 
This oscillation is due to gas in the Hill sphere. For each case, 
the average is plotted as a straight horizontal line. Both the mean 
torque and power vary by 10 percent with varying p. So does the 
corresponding a, while e varies by about 5%. 

A fifth curve, labeled p = 0, is drawn on each panel of 
Fig.|4]which corresponds to the case where / = 1 (no tapering). 
Under these conditions, with an eccentric orbit and an irregu- 
lar, non-stationary density structure (also within the Hill sphere), 
the measured torque and power are quite different, and somehow 
noisy. The amplitude excursions are so large (into the negative) 
that the mean values for the p = case, averaged over one or- 
bit, give a/p = -10 3 instead of 2.8 x 10~ 4 with p = 0.8, and 
tip = -7.14 x 10~ 4 instead of -7.95 x 10~ 5 . This clear differ- 
ence of the p = case, and the very good agreement between 
the other four cases p e [ 0.6 ; 0.9 ], enlightens the importance 
and necessity of tapering. 

For a better understanding of the process, on Fig.|5]are drawn 
in the frame rotating at constant angular velocity fio the traces 
during one orbit of the vectors r p , v^Inr, v^InF, and F (respec- 
tively the position of the planet, its velocity in the non-rotating 
and rotating frames, and the force it feels from the disc), in 
the normalised units. The force is measured with p = 0.8 and 
is magnified by a factor 5 x 10 6 for visibility. The vectors are 
drawn at the beginning of the orbit; then, every twentieth of 
orbit, a cross symbol is drawn on the curve. In the frame ro- 
tating at the velocity Qo> the planet describes an epicycle cen- 
tred on (x - l,y = 0), in the clockwise direction starting at 
(x = a + e — 1 . 15, j = 0). The planet velocity describes also an 
epicycle, centred on (v x = 0, v y - a£2o), in the anti-clockwise 
direction, starting at (v* = 0, v y < 1) ; the velocity in the rotating 
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Fig. 4. Specific torque (left panel) and power (right panel) acting 
on the planet during one orbit for a fixed semi-major axis and 
eccentricity, e = 0.15. The different curves correspond to mea- 
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horizontal lines correspond to the average in each of the four 
cases p + 0. 
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Fig. 5. Evolution during one orbit of the vectors r p : position of 
the planet in the frame rotating at velocity Qo (right) ; v p |nr : 
velocity of the planet in the non-rotating frame (top) ; v p |q F : ve- 
locity of the planet in the frame rotating at constant velocity fio 
(middle) ; F : force from the disc on the planet (left, arbitrary 
unit). 



frame v p \qf is plotted as a dotted line, describing a curve around 
zero in the clockwise direction. As could be expected, the force 
felt by the planet is directed toward the inner disc (F x < 0), more 
precisely in the direction of the wake. The wake tends to rotate at 
a constant velocity. So, between apoastron and periastron, when 
y p is negative (the planet is late with respect to the average ro- 
tation at speed Qo), the wake leads the planet and F y is positive. 
After periastron, the planet leads the root of wake, y p > and 
F y < 0. This can be seen on the density maps of Fig. [6] This 
explains the variations of the sign of the power v^Inr.^ observed 
in Fig. [4] The torque r p A F is always positive because the angle 
(r p ,F) always remains slightly smaller than n, as can be checked 
on the figure. 
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Fig. 6. Density maps, in linear grey scale, in the frame rotating 
at constant velocity Q ( > at four different moments of the orbit 
corresponding to Fig. |4] 

2.3. Discussion 

In this section, we have treated only a particular case, with a 
given planet mass and given disc parameters. The width, the 
depth and the shape of the gap opened by the planet would 
change if the disc viscosity, aspect ratio, or the planet mass vary, 
and the inner disc effect would consequently differ. Also, some 
numerical issues, like the smoothing length for the potential, the 
tapering, and the choice of the frame can affect the measured 
force from the gas on the planet, from a quantitative point of 
view. As seen from the right panel of Fig. [4] a small change in 
the power, e.g. due to numerical aspects, may possibly lead to 
a sign reversal of the average value (implying a reversed migra- 
tion). Exploring all the parameter space would be prohibitive and 
is beyond the scope of this paper. We simply wanted to demon- 
strate here that the inner disc has a non negligible damping effect 
on eccentric giant planets, and we provide a qualitative explana- 
tion of this phenomenon. 

During the evolution towards the final equilibrium state the 
mass of the disc is decreasing. At the end of the simulations it 
lies between 30 and 50 percent of the planet mass in all cases. 
Then the question arises whether one could have some eccen- 
tricity excitation in the disc induced by the planet through a tidal 
instability mechanis m operating th rough the inner 3:1 Lindblad 
resonance in the disc (lLubowll 991). In case of the inverse situa- 
tion, i.e. a planet orbiting inside an outer disc, it has been shown 
that for sufficiently large planet masses (M p Z 3Mj up ) the outer 
disc will turn eccentric even in case of a planet o n a circular or- 
bit jPapaloizo u et al . 200ll lKlev&Dirks en 2006). The strength 
of such a mechanism scales inversly with the planet mass such 
that for the smallest planet masses the timescale of eccentricity 
growth is several thousand planetary orbits. For the instability 
to work the gap created by the planet must be wide enough such 
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name 


M. sin i 


period 


a 


e 




(M Jup ) 


(days) 


(AU) 




b 


1.935 


60.93 ± 0.03 


0.20783 


0.029 ± 0.005 


c 


0.56 


30.38 ± 0.03 


0.13 


0.218 ± 0.002 


d 


0.023 


1.94 


0.020807 


0. 



Tabl e 1. Parameters of the system GJ 876 as given by 
lLaughlin et ail d2005l) for GJ 876 b and GJ 876 c, and by 
http : //exoplanets . eu/ for GJ 876 d. 



that the outer 2: 1 Lindblad resonance, which damps eccentricity, 
has been cleared. Adding a small planet ary eccentricity will r e- 
duce the necessary planet mass slightly (D' Angelo et al. 2006). 

To check for this effect in our simulations we have contin- 
ued models for e = 0.15 for over 3500 orbits and have not seen 
any indication for eccentricity growth in the disk, even in a test 
case where the planet mass has been increased to 5Mj up . We can 
conclude that for our planet masses we do not expect an eccen- 
tric inner disc. Anyway, our simulations have clearly reached 
a steady state when we measure the torque and power of the 
force of the inner disc on the planet, and we expect eccentricity 
damping of the planet by the nearly circular disc. As pointed out 
above, we checked this behaviour through simulations where we 
release the planet from its fixed orbit and follow its orbital evolu- 
tion self consistently. Here we find exactly the predicted results 
for the change in a and e of the planet. 

The non occurence of an eccentric inner disc can be ex- 
plained by the effect that due to the small mass of the planet 
and the not small viscosity and pressure, not all the eccentricity 
damping resonances in the disc are cleared, even for the large 
eccentricities of the planet. 

3. Application to the GJ 876 system 

3. 1 . Characteristics of the system 

GJ 876 is a O.32M M4 star that hosts 3 planets, the two largest 
of them being in a 2:1 Mean Motion Resonance (MMR), see 
TableQ] The resonant angles that characterise this resonance are 
9\ — A c — 2Ab + oj c , 02 = A c - 2Ab + o>b, and Ao> = o>b - cj c , where 
A is the mean longitude and co i s the longitude of the periastron. 
They all three librate around 0°. Laughlin et al. (2005) estimated 
the amplitude of their librations : |0i| max = 7 ± 1.8°, l^lmax = 
34 ± 12°, and |Aw| max = 34 ± 12°. 

Like all the hot Jupiters, the two giant planets in GJ 876 most 
likely formed further out in the protoplanetary disc (beyond the 
snow-line, where water is solid and can contribute to the for- 
mation of a massive core) and then migrated toward the central 
star. Migration is also supported by the resonant motion : such a 
configuration can only be reached by convergent motion of the 
planets. Thus, the outermost planet must have been migrating in- 
ward faster than the innermost one and eventually captured it in 
resonance ; or both planets may have been orbiting in a common 
gap in the disc, repelled one onto the other by the inner and outer 
discs. Then, they can migrate together in this configuration until 
their current position. Dissipation in the gas disc can also modify 
the amplitude of libration of the resonant angles. 

The migration o f two giant planets in MMR i n a gas disc has 
been studied first by Masset & Snellgrove (2001 ) for the case of 
Jupiter and Saturn, and in more details by iMorbidelli & Cridal 
(2007). They showed that if the outermost planet is significantly 
lighter than the innermost one, the migration of both planets may 
well be stopped or reversed outward. In our case of GJ 876 the 
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outer planet is substantially more massive than the inner one, 
so that we expect an inward migration of the pair of planets. 
Under this consideration, the observed small semi-major axes 
of GJ 876 b and GJ 876 c are not a surprise. However, the in- 
nermost of the two resonant planets (GJ 876 c), being pushed 
inward by the outer one (GJ 876 b), should have its eccentricity 
dramatically raised, according to N-body simulations with arti- 
ficial migration rate and detailed hydro-dynamical simulations, 
as discussed in the Introduction. But we have seen that an inner 
gaseous disc has a damping effect on the eccentricity of a giant 
planet, at least for e > 0.1. Thus, we suggest that the presence 
of a gaseous disc inside the orbit of GJ 876 c could prevent its 
eccentricity from reaching unreasonable large values. To verify 
this hypothesis, below we compute numerical simulations of the 
GJ 876 system embedded in a disc, using a hydro-code. 

3.2. Code description and numerical parameters 

Numerical scheme : To have the whole disc simulated, we used 
the cod e bv lCrida et all (120071) . derived from FARGO bv lMassetl 
(2000a,b). A 2D polar grid covers the region where the plan- 
ets orbit, surrounded by a ID grid that extends over all the disc, 
the disc being assumed axisym metrical far from the planets. As 
pointed out in lCrida et alJ d2007l) . this ID grid enables us to take 
into account the global disc evolution, which governs the type II 
migration of the giant planets. In addition, the inner disc evolu- 
tion is self-consistently computed down to an arbitrarily small 
inner radius, that could not be reached by a 2D grid for numeri- 
cal reasons. Thus, the planets can feel the influence of a realistic 
inner disc. Consequently, this code is well adapted to the prob- 
lem that we wish to study. 

In contrast to the previous calculations where we have used a 
non-inertial frame centred on the star, the usage of the added ID- 
grid requires that the frame is inertial and centred on the centre 
of mass of the system (here : star + planet + disc). 

The tapering function used here is / as given in Eq. (0, with 
p = 0.6. 

The 2D grid spans over the region where the planets orbit : 
r e [0.055 ; 0.655], with a resolution of N s - 500 sectors in 
azimuth and N r = 300 rings in radius ; the rings width is thus 
6r = 0.002 AU, which also applies for the ID grid. The outer 
edge of the ID grid is arbitrarily fixed at r — 10 AU, which is far 
enough. The inner edge will be discussed below. 

Disc parameters : The damping effect of the inner disc is pro- 
portional to its mass, in particular to the gas density and the 
amplitude of the wake close to the planet. The deeper the gap 
opened by the innermost planet, the smaller the wake there. The 
wider the gap, the further the disc lies from the planetary or- 
bit. Consequently, a larger gap leads to a smaller damping. The 
shape of the gap is determined by the gas parameters v and h 
dCrida et al.ll2006l) . Thus, these parameters play a crucial role 
for the damping as well. Here, the gas viscosity is still given 
by an a-prescription, with a = 10 2 . The chosen aspect ratio is 
h = 0.07. 

Radiu s of the inner edge of the disc: iCrida & Morbidellil 
(120071) have shown that the inner disc evolution is strongly de- 
pendent on the radius of the inner edge of the disc, and more 
precisely on the ratio between this radius and the radius of the 
planetary orbit. Generally, this radius is poorly constrained, and 
strongly model-dependent. So, the mass of the inner disc could 



be very uncertain : shall this radius be close to 0. 1 3 AU, and there 
would be no inner disc in GJ 876 ; shall it be close to the stellar 
radius and there would be room for a massive inner disc. 

Fortunately, in the case of GJ 876, a 'hot Neptune' is present 
at 0.02 AU. This gives a strong constraint for the location of 
the inner boundary. Indeed, the migration of this planet stopped 
there for some reason. 

A first possibility is that this planet, not massive enough to 
open a gap, an d thus migrating in the type I regime, was caught 
in a planet trap (Ma sset et alj2 006). This planet trap could be the 
inner edge of the disc : indeed at this location, the disc density 
increases rapidly from zero, leading most probably to a positive 
gradient of the vortensity and a strong corotation torque. Then, 
the radius of the inner edge of the disc should have been close to 
0.02 AU. 

A second possibility is that GJ 876 d migrated inwards in 
the gas disc until it lies in the empty cavity inside the edge of 
the disc. A few reasons may explain that this planet crossed the 
planet trap created by the disc edge. The trap may not have been 
strong enough (a jump in density over a only few scale heights 
of the disc is required) ; also, the aspect ratio and the viscosity of 
the disc may have been so low at this place that this planet could 
perturb the disc profile destroying the trap. Finally, turbulence in 
the disc causes random torques, helpful to jump the trap. Once 
the planet is in the cavity, it feel s a negative torque from the 
disc at outer Lindblad resonances (Goldreich & Tremaine 119791 
1980). Thus, the planet goes on migrating inward until there is 
no more gas at the location of its 2:1 resonance (the outermost 
one). In that case, the inner edge of the disc must have been 
located at the 2:1 resonance with GJ 876 d, that is at 0.033 AU. 
Then, as the disc is evaporated by the star from inside out, the 
planet remains there. 

Let us focus on this second possibility. We claim that the 
inner disc could not extend further inward than 0.033 AU, oth- 
erwise it would have pushed GJ 876 d inward, and that it should 
have extended exactly down to this radius, otherwise there would 
be no reason for GJ 876 d to migrate inward to its present posi- 
tion. Thus, the open inner edge of the ID grid will be located at 
r = 0.033 in our simulation. 

3.3. Results 

The innermost planet, GJ 876 d, is not computed. The two largest 
ones are launched on circular orbits at r = 0.36 and 0.21 AU re- 
spectively. At first, the planets influence each other and influence 
the disc, but do not feel the disc effect. So, the planets shape a 
gap in the density distribution and the gas disc reaches an equi- 
librium for this planetary configuration. This phase lasts for 75 
years, which makes about 200 orbits for the outer planet. The 
initial density profile as well as the profile at this time are shown 
on Fig. Q; the final density profile is also plotted. The mass of 
gas present in the 2D grid after this first phase is 2.69 x 10~ 2 
stellar mass (~ 1.7 x 10 28 kg). 

Then, the planets are released and allowed to move under the 
influence of the gas. The evolution of their orbital elements is 
shown on Fig. [8] As expected, the outer planet migrates inward, 
pushed by the outer disc (curve labelled «2 on Fig. [8}. However, 
the inner planet also migrates inward, although less rapidly. This 
is because it did not open a very clean gap (see Fig.|7| and it lies 
on the inner edge of the gap opened by GJ 876 b ; consequently 
the inner planet feels a strong negative corotation torque. 

The 3 relevant resonant angles associated with the 2: 1 reso- 
nance are shown on Fig. [9] After ~ 110 years, the outer planet 
catches the inner one in its 1:2 Mean Motion Resonance (6\, 
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Fig. 7. Density profiles at time (dashed line), at the moment 
where the planets are released (solid line), and at the end of the 
simulation (after 250 years, dotted line). 
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Fig. 8. Semi major axis (blue, left y-axis) and eccentricity (red, 
right y-axis) evolution of GJ 876b (light colour, a\, e\) and 
GJ 876 c (dark, a 2 , ez). The horizontal lines correspond to the 
observed values, as given by Table [TJ 
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Fig. 9. Resonant angles associated with the 2:1 resonance 9\, 
Aco as functions of time. 
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semi major axis evolution of the two planets is almost not af- 
fected ; indeed migration is dominated by the outermost planet, 
which is pushed inward by the outer disc, and that pushes inward 
the inner planet through the resonance locking. Consequently, 
the curves of a,-(f) overlap the ones of Fig. [8] On the contrary, the 
behaviour of the eccentricity of the inner planet e i changes dra- 
matically. It increases faster and continuously, as expected from 
N-body simulations. In fact, both eccentricities rise up to high 
values that are not compatible with the observations {e\ ~ 0.45, 
e2 > 0.1). The eccentricities evolution is displayed on Fig.fTOlfor 
the previous case (labelled "ref"), and in the case where the ac- 
tion of the disc on the inner planet is switched off (labelled "no 
inn.")- This convincingly shows the very important role of the 
inner disc in eccentricity damping for the inner planet, which in 
turn affects the eccentricity of the outer one as well. 



02 and Alo start librating around 0° with small amplitude), and 
the pair of planets goes on migrating in this configuration. The 
eccentricities rise, as expected. But after a phase of eccentric- 
ity growing, a limit value is reached for e\ and ei at about 150 
years. The planets go on migrating at the same rate, while their 
eccentricities remain constant. The value obtained for the eccen - 
tricities is very close to the one given bv lLaughlin etail (120051) . 
After ~ 270 years, the planets have reached their present semi 
major axes, and their eccentricities are oscillating in the ranges 
0.019 <e 2 < 0.032 and 0.21 < e x < 0.25. The amplitude of the 
libration of the resonant angles at this point are \Q\ | max « 18° (but 
it was about 8° at 220 years), |0 2 |max ~ 28°, and |Aw| max * 20°. 
These amplitudes are somewhat larger than the ones given by 
Lau ghlin et al.1 (120051) . but the agreement for the semi majors 
axes and eccentricities is excellent. 



Role of the tapering function : We should also mention that 
in the standard simulation, if one takes p = 0.8 instead of 0.6 
in the tapering function / (see Eq. (O and Fig. the inner 
planet reaches a higher limit eccentricity (between 0.275 and 
0.3), while the eccentricity of the outer one saturates between 
0.035 and 0.05. The eccentricities evolution in this case is also 
displayed on Fig.[T0l(curves labelled "p = 0.8"). The eccentric- 
ities still do not reach extremely high values thanks to the disc 
damping, but this is not in very good agreement with the obser- 
vations (horizontal lines on the figure), especially for what con- 
cerns the inner planet. This also shows that the tapering function 
can have an influence on the final eccentricity of the planets in 
numerical simulations (as could be expected from Fig. and 
one should take this into consideration. 



3.4. Discussion 

Role of the inner disc : For comparison, the same simulation 
has been run with no action of the disc on the inner planet. The 



Disc dispersal : To match with the present state, one has to dis- 
perse the disc when the planets have reached their present semi- 
major axis (at time t ~ 270 yr). This is a common problem when 
modelling the extra-solar planets. 
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Fig. 10. Eccentricity evolution of GJ 876 b (light colour, e<£) and 
GJ 876c (dark, e\) in three cases : (i) standard simulation (same 
as Fig. [8] curves labelled "ref") ; (ii) same as (?) but the inner 
planet is not affected by the disc (curves labelled "no inn.") ; (Hi) 
same as (/) but the p parameter of the tapering function / is 0.8 
(curves labelled "p = 0.8"). The horizontal lines correspond to 
the observed values, as given by Table Q] 



To check what happens when one removes the gas dis c, we 
applied a procedure a lso used in Morbid eTli et alj d2007l) and 
Thom mes et al] (|2007) to have the gas disappear smoothly for 
the planetary system not to be perturbed by a sudden change in 
the potential : from time t = 250yr, the gas density is decreased 
exponentially in each cell of the grid with a time scale of 27.5 
years. The result is shown on Fig. QT|: as the gas disappears, 
the planets remain in 2: 1 MMR, while their migration speed de- 
creases exponentially (it is well known that when the planet is 
more massive than the disc, the inertia of the planet is the lim- 
iting factor in the type II migration regime). One should note in 
passing that as the migration speed and the eccentricity damping 
are both proportional to the gas density, the K factor is not af- 
fected in this procedure ; consequently, the equilibrium value of 
the eccentricities is not affected, and they remain close to 0.03 
and 0.22 respectively, while the semi major axes converge to 
0.21 and 0.13 AU. In the end, gas has almost completely dis- 
appeared, and we are left with a planetary system very similar to 
GJ 876. 

The disc clearing process is complex and not well con- 
strained, in particular in the vicinity of the star. In general, the 
gas density first slowly decreases while the disc accretes onto the 
central star and spreads outward. When the density is low, the 
photo-evaporation by the central star can play a significant role 
and erode the disc. The extreme UV photons ionise and heat 
the upper layer of the disc, so that gas can leave the potential 
well of the star, and the density decreases. However, this works 
more efficiently at radius larger than about 1 AU (close to the 
star, the gravity is too strong). Consequently, the region where 
the two giant planets of GJ 876 are orbiting should not be much 
affected ; in particular, the inner disc should not disappear. The 
remnant disc inside ~ 1 AU will viscously spread onto the star 
and out. So, the migration path of the planets should not be much 
affected. In addition, the viscous evolution dominates the disc 
evolution until the density is very low. Then, photo-evaporation 
happens on a timescale that is much shorter than the disc life 
time. Thus, the final phase of gas dispersal occurs when the disc 
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Fig. 11. Semi major axis and eccentricity evolution of GJ 876 b 
and GJ 876 c, with disc clearing from time t - 250yr on. The 
horizontal lines correspond to the observed values, as given by 
TableG] 



is too low to have any significant influence on the planets on such 
a short timescale. 

Thus, we think that the planetary configuration should not 
be significantly perturbed during this phase, and that the above 
modelling by exponential damping of the density gives reliable 
results. Anyway, these planets somehow migrated toward their 
host star, until the gas density was too low to push them. They 
were left at 0. 13 and 0.21 AU, and then the disc disappeared. Our 
simulation shows that when the planets stop migrating, they au- 
tomatically have the correct eccentricities ; hence, our simulation 
is consistent with the observed configuration. To our knowledge, 
this is the first time that an extra-solar system is reproduced in a 
fully hydro simulation taking into account all the protoplanetary 
disc and allowing for a significant migration of the planets with 
correct eccentricities. 



4. Modelling an inner disc by N-body calculations 

Studying the effect of an inner disc on the evolution of the reso- 
nant planetary systems by full hydro-dynamical calculations (as 
done in the previous section) requires typically a large amount 
of computer time. Fortunately, the effects of the outer and inner 
discs can also be modelled approximately conveniently by grav- 
itational N-body simulations using properly parameterised non- 
conservative forces. These forces can be derived by using the mi- 
gration rate a/ a and the eccentricity dampin g rate e/e of a planet 
embedded into the protoplanetary disc (see Lee & Pe ale 2002a; 
Beauge et al. 2006, for two different approaches). Instead of the 
migration and damping rates one can also use the corresponding 
e-folding times defined as r a - -(a/ay 1 and j e = — (e/e)~ l . 

When studying the formation of a resonant system consist- 
ing of an inner and an outer giant planet, usually the outer planet 
is forced to migrate inward. When the ratio of their semi-major 
axes approaches a critical limit, a resonant capture may take 
place between the m (for the conditions o f a resonant capture into 
the 2: 1 MMR see lKlev & Sandorll2007l) . After the resonant cap- 
ture the two planets migrate inward as the outer one still feels 
the negative tidal torques of the disc. 

Having studied in the previous section the hydro-dynamical 
evolution of GJ 876 thoroughly, in this section we provide fur- 
ther results obtained in the framework of the three-body problem 
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with dissipative forces. Our results show that the presence of an 
inner disc during the migration of the giant planets is consistent 
with the observed state of the resonant systems having two giant 
planets engaged in a 2:1 MMR. 

Since the eccentricity of the inner planet is excited by the 
resonant interaction, its orbit becomes more and more elongated 
penetrating into the inner disc. This represents a damping mech- 
anism which acts against the eccentricity excitation, and may set 
a quasi-equilibrium between these processes keeping the eccen- 
tricity of the inner planet on a constant value (within certain lim- 
its) during the whole migration process. The effect of an inner 
disc can be investigated by using a repelling non-conservative 
force acting on the inner planet parameterised by a positive mi- 
gration rate a/a and a negative e/e. Similarly to the case of the 
inward migrating outer planet, a ratio K of the above parameters 
can also be defined. According to the definition of the e-folding 
times, t„ will have a negative sign in this case. 

We note, however, that when both the outward and inward 
migration of the outer and inner planets are considered simul- 
taneously, the final state of the system cannot be characterised 
uniquely by using only the K ratios for the inward and outward 
migration. The final values of the semi-major axes and the ec- 
centricities depend directly on the migration rates (or on the cor- 
responding e-folding times) d\/a\, 02/02 and the eccentricity 
damping rates e\je\, <?2/<?2 (where indices '1' and '2' stay for 
the inner and outer planet, respectively), and not only on their 
ratios K\ and A^. The characterisation of the system's final state 
by using the migration parameters will be discussed more de- 
tailed at the end of this section. 

In the f ollowing we r epeat t he three-body calculations for 
GJ 876 by iLee & Peald d2002al) adding the effects of an in- 
ner disc, then we review the simulations for HD 73 526 by 
Sand or et alj (|2007), and finally present our new results in mod- 
elling the formation of HD 82 943 and HD 128 31 1 with an inner 
disc. 

4.1. GJ876andHD 73526 

GJ 87 6 : First we repeated the simulations of ILee & Peale 
(2002a) on the formation of GJ 876 by using three-body cal- 
culations with dissipative forces. Similarly to their results, we 
required K2 = T a ,/T e2 = 100 if only the outer giant was affected 
by an outer disc. The evolution of the semi-major axes and the 
eccentricities is shown in the top panel of Fig. [12] In our calcula- 
tions we used r 02 = 2x 10 4 and T e2 — 2x 10 2 years giving exactly 
K2 = 100. We note again that this ratio between the e-folding 
times is very high and maybe physically unrealistic for massive 
planets. 

To quantify the damping effect of the inner disc on the inner 
planet while the giant planets revolve in a common gap, we per- 
formed further three-body simulations with also r fl] = -2 x 10 4 
and T e , = 2.5 x 10 3 years, where we point out that the minus sign 
of T fll stands for the outward migration. To model the damping 
effect of the outer disc on the outer planet we used the e-folding 
times T flJ = 2 x 10 4 and r f , = 2.5 x 10 3 years. Note that the 
above migration parameters correspond to K\ = K2 = 8. The re- 
sult of our calculations is shown in the bottom of Fig. [12] These 
figures in fact are very similar to those obtained bv lLee &~ Peale 
d2002al) . however to model the behaviour of the system, we used 
different migration parameters. 

HD 73 526: As we already menti oned, the effe c t of a n inner 
disc was studied more detailed by ISandor et alj d2007l) in the 
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Fig. 12. Behaviour of the semi-major axes and the eccentricities 
of the resonant giant planets in the system GJ 876. The horizon- 
tal lines correspond to the observed values of the eccentricities. 
Top: Only the outer disc is taken into account, and therefore 
only the orbital evolution of the outer planet is affected, with 
the e-folding times r fl2 = 2 x 10 4 years, t Ci — 2 x 10 2 years, 
so K2 = 100. Bottom: Same as above in the presence of outer 
and inner discs. In this case, T ai = — 2 X 10 4 , T e , = 2.5 x 10 3 , 
T ai = 2 x 10 4 , and r e2 = 2.5 x 10 3 years (giving K X =K 2 = 8). 



case of HD 73 526. There, it has been shown that if the outer 
planet is damped only, one needs K2 = 15 - 20, which accord- 
ing to the hydro-dynamical simulations seems to be a too high 
value. Moreover, if only the outer planet is damped, the eccen- 
tricity of the inner planet shows an increasing tendency, which 
can result in exceeding the limit coming from the observations 
(e; nn as 0.3). On the other hand, if in addition the inner planet is 
damped by an inner disc, the eccentricities will level off and the 
outcome of the migration scenario fits better to the observations : 
the formation of HD 73 526 could be modelled successfully by 
using the e-folding times t Bi = -5 X 10 4 , t«, = 10 4 years and 

= 5 x 10 3 , T ei = 10 3 (corresponding to Ki = K 2 = 10). The 
behaviour of the eccentricitie s in this case is sh own in the top 
panel of Fig. 8 in the paper of Sa ndor et al.l (120071) . 

We can conclude that the presence and effect of an inner disc 
in modelling the formation processes of GJ 876 and HD 73 526 
is physically more realistic, giving results fitting very well to the 
radial velocity observations. In what follows, we show that in the 
case of the systems HD 82 943 and HD 128 31 1 the presence of 
an inner disc is still a real alternative, however these systems do 
not need such strong dampings for the eccentricity of their inner 
giant planets. 



10 



Crida et al. : Inner disc and resonant planets 




10000 20000 30000 40000 50000 

time [years] 




L — 1 1 1 1 1 

10000 20000 30000 40000 50000 

time [years] 



Fig. 13. Behaviour of the semi-major axes and the eccentrici- 
ties of the giant planets in the resonant system HD 82943. The 
horizontal lines correspond to the observed values of the eccen- 
tricities. Top : Only the motion of the outer planet is affected, 
with the e-folding times r fl2 = 2 x 10 4 and T e2 = 2.5 x 10 3 
years (K 2 = 8). Bottom: An inner disc is also assumed, with 
T fll = -2 x 10 5 , r ei — 5 x 10 4 years e-folding times for the inner 
planet (K\ = 4), which allows the use of larger T e , = 2.5 x 10 3 
years and thus lower ratio K 2 = 4 for the outer planet. 

4.2. HD 82943 and HD 12831 1 

HD 82 943 : Recently Eee et all d2006l) presented four sets of 
orbital solutions based on a best-fit double-Keplerian model for 
HD 82 943. Using their orbital elements as initial conditions, di- 
rect numerical integrations show that three of them exhibit or- 
dered behaviour, while one (Fit I in the cited paper) is desta- 
bilised after a few thousand years. In one of the three dynam- 
ically stable orbital solutions (Fit II) the variations of the ec- 
centricities are negligible, and the resonance variables oscillate 
with small amplitudes indicating clearly that the system is deep 
in the 2:1 MMR. Since this behaviour can be reached during an 
inward convergent migration of the giant planets, we investigate 
the formation of this system using the scenario of the planetary 
migration. 

First we assume that only an outer disc is present and the 
outer planet feels its damping effect. To obtain the dynami- 
cal state calculated based on Fit II (ei = 0.422, e 2 = 0.14, 
«i = 0.74 AU, a 2 = 1.18 AU), we needed the ratio K 2 = 8, using 
T ai — 2 x 10 4 and T ei = 2.5 x 10 3 years. This ratio K 2 does not 
seem too high (as it lies typically between 1 and 10), however, 
similarly to HD 73 526, the eccentricities are slightly increas- 
ing during the whole migration process. This is not a problem 
if the migration of the planets is terminated gradually around 
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Fig. 14. Behaviour of the semi-major axes and the eccentricities 
of the resonant giant planets in the resonant system HD 128 3 1 1 . 
The horizontal lines correspond to the observed values of the ec- 
centricities. Top : Only the motion of the outer planet is damped 
by an outer disc, with K 2 - T a2 /T ei - 7. Bottom : An inner disc 
is also assumed. The e-folding times for the inner planet are 
T a , = 2 x 10 4 , T e , = 10 4 years (Ki = 2), which allows to use 
much lower ratio K 2 -2 for the outer planet. 



the observed values of their semi-major axes but it may result 
in exceeding the limits derived from observations, see the top 
panel of Fig. [13] On the other hand, by assuming the presence 
of an inner disc, the eccentricities reach their constant values 
gradually which do not seem to change further during the migra- 
tion process. The latter behaviour is shown in the bottom panel 
of Fig. [13] During this simulation we used for the inner planet 
r a , = -2 x 10 5 and r e , = 5 x 10 4 years, and for the outer planet 
T m - 2 x 10 4 and T e2 = 5 x 10 3 years (giving Ki - K 2 - 4). Since 
the eccentricities seem to reach constant values, this case seems 
to be a more convenient formation scenario for HD 82 943 than 
the previous one. 

We can conclude that the present dynamical behaviour of th e 
resonant system HD 82 943 (based on Fit II of lLee et alj d2006h ) 
can be explained in both ways : by assuming an outer disc alone 
or by assuming an inner and an outer disc. However, the latter 
case seems to be more reasonable since when the migration oc- 
curs over longer times it guarantees constant eccentricity values 
during the whole process. 

HD 128311 : Finally, we present our results for HD 128 311. 
The most recent orbital solution for this system is presented by 
IVogt et all d2005l), while poss ible formation scenarios are out- 
lined bv lSandor & Klevl d2006l) . In the latter study, the formation 
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of the resonant system is modelled by an inward convergent mi- 
gration of the giant planets which is followed by a sudden pertur- 
bation ; only the effect of an outer disc was considered. Prior to 
the sudden perturbation, the outer giant planet migrated inward 
and captured the inner one into a 2:1 resonance. It is essential 
to find the final values of the eccentricities as due to the mi- 
gration process since the perturbative events modify only their 
oscillations. These values are e\ m 0.45 and e2 ~ 0.15 being 
obtained after a migration characterised by Ki = T fl ,/r e , = 5. In 
the present study we need somew hat larger ratio, K2 = 7, since 
contrary to lSandor & Klevl d2006). we do not stop the migration 
when the planets reach their actual positions. The evolution of 
the semi-major axes and eccentricities are presented in the top 
panel of Fig. [14] 

Assuming an inner disc, the same final state of the migratory 
evolution of HD 128 31 1 can be obtained. In our numerical ex- 
periments we used r fl] = -2 x 10 4 , r e , = 10 4 years for the inner 
planet, and r fl2 = 2 x 10 4 , r ei — 10 4 years for the outer planet 
(being equivalent to K\ — K2 = 2). The evolution of the system 
is shown in the bottom panel of Fig. [14] Comparing the top and 
the bottom panels of Fig. [14] we can conclude that the formation 
of HD 128 3 1 1 can also be modelled by assuming an inner disc, 
however in this particular case, this is not definitely necessary. 




Fig. 15. Pairs of migration parameters t , and r e , which result 
in the observed state of the system GJ 876 for fixed r m = 4 x 
10 4 years and r e2 — 5 x 10 3 years. Crosses : Empirically found 
with N-body simulations. Dashed curve : Analytical expression 
Eq. ( [TBI (see Sect. 14.4) . Solid line : Constant ratio K\ = 8, for 
comparison. 



4.3. Final state of a migration with an inner disc: Numerical 
experiments 

In this part we show the results of additional numerical simula- 
tions in which we studied how the different parameters of migra- 
tion influence the final state of the resonant system. It was quite 
evident from the beginning that the characterisation based only 
on the ratios K t and K2 would not yield unique results in the case 
when the inner planet is also affected by an inner disc. 

Our calculations are based on the observed orbi t al and p hys- 
ical parameters of GJ 876, as given bv lLee & Pea le (2002a). For 
the inward migration of the outer planet we fixed r fl , = 4 x 10 4 , 
T ei = 5 x 10 3 years (so, K2 - 8), and we changed t 0j and T e , in 
such a way to obtain the observed state of the system. 

As starting values to our numerical simulations we used 
r a , = — 3 X 10 4 years and found that r e , = 3.7 x 10 3 years gives 
a correct result (Ki » 8 here). Then we increased gradually the 
absolute value of r fl] corresponding to a weaker damping on a\. 
To obtain the observed behaviour of the system, for larger r a , , 
we needed larger r ei meaning a weaker damping on e\ too. In 
order to explore the mutual dependence of T a , and T e] , we per- 
formed a series of numerical experiments. Our results are shown 
in Fig.1151: on the jc-axis the absolute value of r 0l , while on the 
y-axis the eccentricity damping time r e , are displayed in loga- 
rithmic scale. The crosses show those corresponding values of 
T fl , and T e , which are needed to obtain the same final values of 
the eccentricities e\ and <?2 (0.26 and 0.035 respectively). The r e , 
increases rapidly for small |r fl| |, however it is not proportional to 
r a , (compare with the straight line on Fig. [TBI ), but it levels off 
and clearly tends to a limit r fimax , which in this case is slightly 
higher than 9 x 10 3 years. 

We can conclude that for a fixed pair of r a , and T ei , there is 
no unique K\ that determines the final state of the system. On the 
contrary, there exists a r e , max , which determines the final state of 
the system if the inner disc does not affect the semi-major axis 
of the inner planet. In reality, however, the inner disc (as rotating 
faster) can transmit angular momentum to the inner planet as 
well as energy, increasing its semi-major axis. In this case we 
need smaller r ei , or equivalently, more effective damping on e\. 



If we compare the two migration scenarios : (i) only an outer 
disc is present and there is no inner disc between the inner planet 
and the star ; and (ii) both an outer disc and an inner disc are 
present and the planets orbit in a common gap between them, we 
can summarise the following results : The final state of the sys- 
tem in case (i) can be described solely by the ratio K2 = T a2 /r e2 . 
In case (ii), the ratios K\ and K2 are not enough, the final state 
of the system depends on additional parameters too ; in the most 
general case, when the inner disc forces the inner planet to mi- 
grate outward, there are essentially three parameters, which can 
be K\, K2, and the ratio r fl] /r a , for instance. In next subsection, 
we analyse theoretically this behaviour. 



4.4. Final state of a migration with an inner disc : Analytics 

For the semi major axis and eccentricity of a planet to de- 
crease exponentially with damping timescales t„ = -dja and 
T e = —e/e, its energy E and angular momentum H must vary 
with the following rates (through Eqs. © and 0) ) : 



E/t c1 
H 

H = — . 
2 \1 



2e 2 1 



(6) 
(7) 



When two planets are considered, the variation of the energy 
of the system { pair of planets }, E is the sum of the energy vari- 
ations applied to both planets, so that 

. Ei Ej 

£ = — + — . (8) 

If the two planets are in resonance, the ratio between their semi 
major axes is kept constant, and consequently also the ratio be- 
tween their energies. Let us define : 



s = E 2 /Ei = M 2 ai/Mia 2 . 
Then, 

E = E x +E 2 = (1 +s)Ei , 
E = (1 +e)Ei . 



(9) 

(10) 
(ID 
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Note that due to some energy exchange between the planets 
through the resonance, the variation of the energy of the first 
planet is not Ei/r ai , but E\ — E/(l + s), with E given by Eq. jSJ. 

The same holds for the angular momentum. For the angular 
momentum H of the system { pair of planets } (H - H\ + H 2 ), 
the total variation rate reads (from Eq. (01 ) : 



^,211-^ tJ 



(12) 



For the final state, when the planets are in resonance and their 
eccentricities are also constant, as is the ratio of their angular 
momentum. Let us define : 



„ ,„ M 2 a 2 (l-e 2 2 ) 
77 = H 2 /Hi = — 



M[ \ai(l-ei 2 ) ' 
Then, H = H\ + H 2 = (1 + tj)H u and also 
H = (l+t])Hi ■ 



(13) 



(14) 



If the eccentricities are constant (e, = 0), then, from Eqs. (O, 
COB, CED, and finally © : 



' ~ 2 E, ~ 2 E ~ 2 I 1 + e 



(15) 



Using Eqs. ( TTZb and (IT~5b . Eq. ( TPfl i gives a relation between 
T fll , Ta 2 , r ei , and t C2 . In the problem that we study, the masses of 
the planets are known, as well as their eccentricities and the ratio 
between their semi major axes. The question is : for a given effect 
on the outer planet, what should the effect on the inner planet 
be to be consistent with the observed eccentricities ? Below, we 
solve this equation in the unknown r e , with the free parameter 
T fll , while T fl2 , T e% , s and 77, are given constants (e and 77 being 
defined by Eqs. © and ( [T3l respectively). Using Eqs. ( fT2b and 
( TBI Eq. (Hi reads : 



Hi /l/r fll +e/r a 

-(l +'7)-r — n — 

2 \ 1 + e 



Zf 

/=1,2 Z 


1 2e? 1 




2ej 2 1 




1-ei 2 r ei 


- — +77) 


2 ei 2 1 


1 f ( 


1-ei 2 r ei 


in 

Ta 2 \ 



+ — 1- 



Eventually, we obtain : 
1 l-e\ 2 [ 1 77 - e 1 2e 2 2 rj \ 1 1— ei 2 e — 77 



2ei 2 \T a , 1 + e r C2 1 - e 2 2 ) T ai 2e\ 2 1 + e 



(16) 



1/T„ 



As one can see, the damping rate that should be applied to 
the inner planet (1 /r ei ) is the sum of two terms. 

The first one (1/T eimax ) is required to balance the action on 
the outer planet ; indeed, when no force is applied to the inner 
planet, the energy loss rate of the outer planet is not the expected 
E 2 /T ai , but (E 2 /T a2 )-^ because the inner planet also has to lose 
energy to preserve the resonant motion ; thus, the angular mo- 
mentum loss rate is overestimated by the expression Eq. (|7), and 
both eccentricities rise. Therefore, a damping of the eccentricity 
needs to be applied on the inner planet as well. 



The second term (A^o/t^) is proportional to r ai _1 ; the co- 
efficient K\x) is negative if e 2 > (1 - (a 2 la\y'~) + (a2/a\) 3 e 2 2 , 
which is always true in the case of a 2:1 MMR for e 2 < 0.866. 

In the case stud ied in the previous Sect. 14.31 one has 
dLee & Peald l2002al) : M 2 = 1.808M Jup = 5.65 x 10~ 3 M, , 
Mi = 0.5696M Jup = 1.78 x 1(T 3 M« , e x * 0.265 , e 2 * 0.035 
, T fl , = 4 x 10 4 , T ei = 5 x 10 3 and the two planets are in 2:1 
MMR, with a 2 /ai = 1.602. Thus, 77 = 4.1641 and s = 1.9812. 
This gives T e]lnax = 9289 years, and K\fi = -4.8471. The dashed 
curve on Fig. [TBI shows r fl as a function of r fll from Eq. (I161 l. 
The fit is excellent. 

This shows that Eq. (fT6] l efficiently gives what parameters 
one should choose to reproduce a given system by the means of 
N-body simulations with dissipative forces, or one should have 
in the protoplanetary disc. In particular, the expression of r eimax 
shows that if e 2 is relatively small, a huge K 2 of the order of e 2 ~ 2 
is required for 7> imax — > 00. If the inner disc tends to make the 
inner planet migrate outward with r fll < 0, an additional damp- 
ing on its eccentricity is required, which can be expressed as 
T Cl add = T a JKifl, with Ri o < 0, depending on all the parameters 
of the system. On the other hand, if the inner disc tends to attract 
the planet sufficiently (which may happen as has been shown 
in Sect. |2), no eccentricity damping on the inner planet may be 
required at all (this happens in our case if T a , = +45025 years). 

With this equation, one can try to find reasonable parameters 
(with not too big K\ and K 2 ) to explain how all the known reso- 
nant exoplanetary systems can have been formed in the disc. Our 
numerical simulations have shown that this can be achieved for 
at least 4 such systems. 



5. Conclusion 

In this paper, we address the problem of the moderate eccentric- 
ity of extra-solar planets in resonance. The resonant configura- 
tion requires a convergent migration of the planets, but continued 
migration in resonance leads to unlimited eccentricity growth if 
no eccentricity damping mechanism is at work. 

In Sect. [2] we have shown that an inner disc has a non negli- 
gible influence on a giant planet on an eccentric orbit and mod- 
ifies its orbital parameter. The strength of this effect varies with 
planetary eccentricity. For very small eccentricities, e < 0.05, 
we find a very small but positive e while for larger e it becomes 
more and more negative. The induced change in semi-major axis 
remains positive, as expected for the Lindblad torques induced 
by an inner disc on a massive planet. Only for larger eccentric- 
ities e > 0.35 the change in semi-major axis becomes negative, 
leading to inward migration. The usage of a constant /f-factor 
between the e-folding time scales of the semi-major axis and the 
eccentricity is clearly an over simplification. 

However, the measured influence of the disc on the planet 
depends on the adopted tapering function applied to exclude at 
least parts of the gas within the Hill sphere of the planet, and pos- 
sibly on other numerical effects, such as the smoothing length, 
the resolution, or boundary conditions. Hence, determination of 
the exact, absolute magnitude of the effect is a bit difficult. 

Anyway, an eccentricity damping should be applied to the 
inner planet to obtain realistic results, in particular if its eccen- 
tricity is larger than 0.1. Using a hybrid 2D- ID hydro-code that 
allows the simulation of the whole disc from its physical inner 
boundary to an arbitrary large radius, we compute the longterm 
evolution of the GJ 876 system, for which the orbital elements 
are well known, and the radius of the inner edge of the disc can 
be estimated thanks to the presence of a third planet very close 
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to the star. We find that the inner disc does not disappear after 
the planets opened a common gap, and that it effectively helps 
damping the eccentricity of the inner planet. With realistic disc 
parameters (viscosity and aspect ratio), we are able to reproduce 
the observations. 

As hydro-simulations are very time-consuming, we finally 
perform customised N-body simulations with explicit non- 
conservative forces, added to mimic the effect of the outer and 
inner disc. We find that applying a significant damping to both 
planets, as required by the influence of both the inner and outer 
discs, enables us to fit quite well a few other exoplanetary sys- 
tems. In addition, N-body simulations have shown that when two 
planets are considered, the ratio between the eccentricity and 
semi major axis damping applied on each of them (K\ and K2) 
are not sufficient to determine the final state of the system. In 
fact, for a given damping in a and e applied on the outer planet, 
one can express analytically the eccentricity damping that should 
be applied on the inner planet to match the observed orbital con- 
figuration, as a function of its semi major axis damping. 

Given the satisfying fits of a few systems that we obtain, we 
claim that the problem of the low eccentricities of the resonant 
exoplanetary systems simply stems from the fact that the inner 
disc had not been taken properly into account, which is not rea- 
sonable. A pair of planets orbiting in a protoplanetary disc may 
well orbit in a common gap and enter a Mean Motion Resonance, 
but it does not necessarily open a complete cavity from the star to 
the outermost pla net, so that the inner disc s hould influence the 
planets dynamics. ICrida & Morbidelli (2007) already suggested 
that the opening of such a cavity by a giant planet requires very 
specific conditions ; in fact, according to our hydro and N-body 
simulations, the low observed eccentricities of the exoplanets in 
resonance support the idea that the inner disc does not disappear 
in general. 
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